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CONVERSION  FACTORS,  US  CUSTOMARY  TO  METRIC  (SI) 
UNITS  OF  MEASUREMENT 


US  customary  units  of  measurement  used  in  this  report  can  be  converted 
to  metric  (SI)  units  as  follows: 


Multinl 


To  Obtain 


feet 

0.3048 

metres 

feet  per  second 

0.3048 

metres  per  second 

feet-feet  per  second 

0.0929 

metres-metres  per  second 

feet  per  second  per 
second 

0.3048 

metres  per  second  per 
second 

pounds  (mass)  per  square 
foot 

4.882428 

kilograms  per  square  metre 

pounds  per  foot 

1.488189 

kilograms  per  metre 

pounds-second-second  per 
foot  per  foot  per  foot 
per  foot 

52.5540137 

kilogram-second-second  per 
metre  per  metre  per 
metre  per  metre 

EROSION  CONTROL  OF  SCOUR  DURING  CONSTRUCTION 
CURRENT- -A  WAVE- INDUCED  CURRENT  MODEL 


PART  I :  INTRODUCTION 


Statement  of  the  Problem 


1.  Often,  large-scale  engineering  structures  such  as  jetties  and  break¬ 
waters  are  constructed  in  the  nearshore  region  to  stabilize  navigation  chan¬ 
nels  or  protect  harbor  entrances  and  beaches.  The  structures  are  usually 
built  of  rock  or  precast  concrete.  During  construction,  the  massive  struc¬ 
tures  alter  the  waves  and  currents  near  their  location.  Waves  break  on  the 
toe  of  the  structure  and  the  resulting  turbulence  causes  material  from  the 
bottom  to  be  suspended.  The  suspended  material  in  turn  is  moved  away  from  the 
region  by  wave-induced  and  other  currents.  This  results  in  erosion  developing 
at  the  toe  of  the  structure  since  the  lost  material  may  not  be  replaced  by 
natural  processes.  In  order  to  ensure  that  the  structure  will  be  stable  and 
perform  its  function  as  desired,  it  will  be  necessary  to  fill  with  nonerodible 
material  any  scour  hole  that  may  have  developed  due  to  erosion.  As  a  result, 
extra  quantities  of  material  may  be  required  and  construction  costs  may  be 
overrun.  To  minimize  cost  increases  due  to  scour  during  construction,  it  is 
necessary  to  estimate  beforehand  the  likelihood  and  amount  of  potential  scour 
during  construction.  This  is  a  very  complicated  problem  and  the  solution  to 
it  will  depend  strongly  on  the  particular  field  site  and  the  coastal 
environment. 

2.  Since  breaking  waves  and  the  currents  induced  by  them  play  a  vital 
role  in  transporting  sediment  away  from  coastal  structures  resulting  in  scour, 
it  is  important  to  predict  wave-induced  currents,  with  and  without  structures. 
In  view  of  the  environment  and  site-specific  nature  of  the  problem,  a  general¬ 
ized  numerical  model  is  needed.  Results  from  such  a  model  can  be  used  as  input 
to  a  sediment  transport  model  to  predict  erosion  during  construction. 

Purpose  of  the  Study 


3.  The  purpose  of  this  study  is  to  develop  a  generalized  numerical  model 


that  can  predict  currents  induced  by  breaking  waves  at  locations  with  or  with¬ 
out  coastal  structures.  The  model  should  be  applicable  to  real-life  bathym- 
etries  that  are  often  arbitrary  and  irregular  and  should  be  computationally 
efficient  and  economical  in  view  of  the  large  numerical  grids  often  required 
in  engineering  projects. 


PART  II:  THEORETICAL  BACKGROUND 


Review 

4.  In  recent  years,  there  has  been  a  growing  interest  in  the  numerical 
modeling  of  longshore  currents  and  nearshore  circulation  caused  by  the  action 
of  breaking  waves.  Results  of  such  simulation,  besides  being  useful  on  their 
own,  form  an  essential  input  to  numerical  and  physical  models  for  nearshore 
processes  such  as  sediment  transport.  The  design  and  construction  of  large- 
scale  engineering  structures  such  as  jetties  require  the  determination  of  wave- 
induced  currents  not  only  near  the  open  coast  but  also  near  inlets;  in  addi¬ 
tion,  the  effects  jf  the  proposed  or  existing  structures  on  these  currents  are 
required.  Several  publications  on  longshore  currents  and  nearshore  circulation 
have  appeared  in  the  literature  during  the  past  two  decades.  Some  of  these 
relate  to  experimental  or  field  studies  and  others  to  analytic  solutions  and 
numerical  models.  Among  the  latter,  mention  must  be  made  of  Bowen  (1969), 
Longuet-Higgins  (1970),  Thornton  (1970),  Noda  (1974),  Jonsson,  Skovgaard,  and 
Jacobsen  (1974),  Birkemeier  and  Dalrymple  (1975),  Liu  and  Mei  (1976),  Liu  and 
Lennon  (1978),  Liu  and  Dalrymple  (1978),  Ebersole  and  Dalrymple  (1980),  and 
Vreugdenhil  (1980).  Several  of  these  studies  and  models  either  consider  only 
simple  and  idealized  situations,  for  example,  plane  beaches  and  periodic  ba- 
thymetries,  or  neglect  terms  of  the  governing  equations  involving  unsteadiness, 
advection,  and/or  lateral  mixing.  Often,  a  linear  friction  is  assumed.  For 
practical  engineering  application,  a  generalized  numerical  model  is  needed 
that  can  handle  real-life  bathymetries .  The  model  should  be  flexible  in  that 
one  should  be  able  to  change  easily  the  formulation  of  terms  such  as  mixing 
and  friction  as  improved  understanding  of  these  processes  is  gained  in  the 
future.  It  should  be  computationally  efficient  and  economical  for  large 
numerical  grids.  This  report  describes  one  such  numerical  model  for  wave- 
induced  currents  developed  at  the  US  Army  Engineer  Waterways  Experiment  Sta¬ 
tion  (WES).  The  model  called  "CURRENT"  has  been  applied  successfully  to  the 
determination  of  longshore  currents  near  open  coasts  and  wave-induced  currents 
near  inlets.  It  can  handle  impermeable  nonovertopping  structures. 

5.  In  terms  of  its  analytic  formulation  the  present  model  uses,  to  a 
large  extent,  the  approach  of  Ebersole  (1980),  and  Ebersole  and  Dalrymple 
(1980).  Whereas  their  numerical  model  used  an  explicit  method  of  solution. 


the  present  model  uses  a  highly  efficient  alternating  direction,  implicit, 
finite  difference  scheme.  In  view  of  the  similarity  between  the  equations  for 
wave-induced  currents  and  long  waves,  the  present  model  was  created  by  modi¬ 
fying  an  existing,  well-tested  WES  shallow-water  wave  numerical  model  known  as 
WIFM  (WES  Implicit  Flooding  Model)  (Butler  1980) .  The  velocity  version  of 
WIFM  used  here  has  nonlinear  advective  terms.  The  friction  and  mixing  terms 
used  in  WIFM  were  modified  to  conform  to  the  formulations  normally  used  for 
wave-induced  currents.  Radiation  stress  terms  were  added  to  the  model,  since 
these  are  usually  the  "driving"  terms  for  wave-induced  currents.  Certain 
capabilities  of  WIFM  such  as  flooding/drying  and  wind-induced  current  calcula¬ 
tion  were  not  utilized  in  the  model  for  wave-induced  currents. 


Equations  of  Motion 


6.  The  hydrodynamic  equations  used  in  the  model  for  wave-induced  cur¬ 
rents  may  be  derived  from  the  Navier-Stokes  equations  (for  details,  see 
Phillips  (1969)  and  Ebersole  (1980)).  It  is  assumed  in  the  derivation  that  the 
fluid  is  homogeneous  and  incompressible,  and  the  vertical  accelerations  are 
negligible  so  that  the  pressure  distribution  is  hydrostatic.  By  vertically 
integrating  the  three-dimensional  form  of  the  equations  and  applying  appropri¬ 
ate  boundary  conditions,  the  depth-averaged  two-dimensional  form  of  the  equa¬ 
tions  of  motion  and  continuity  are  obtained.  These  equations  are  derived  by 
time-averaging  over  a  time  interval  corresponding  to  the  period  of  the  waves. 
Refc.*.ing  to  a  Cartesian  coordinate  scheme  (Figure  1),  these  are: 
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a.  CROSS  SECTION  A-A 


OCEAN 

BOUNDARY 

UNE 


b.  PLAN 


Figure  1.  Definition  sketch  for  an  irregular  beach 


U  and  V  =  depth-averaged  horizontal  velocity  components  at 

time  t  in  the  x-  and  y-directions ,  respectively, 
ft/sec* 

r|  =  displacement  of  the  mean  free  surface  with  respect 
to  the  still-water  level,  ft 

2 

g  =  acceleration  due  to  gravity,  ft/sec 

2  4 

p  =  mass  density  of  seawater,  lb-sec  /ft 
d  =  r)  -  h  =  total  water  depth,  ft 

h  =  bed  elevation  with  still-water  level  taken  as  zero 
(note  h  is  negative  for  water  cells  and  positive 
for  land  cells),  ft 

t,  and  t,  =  bottom  friction  stresses  in  the  x-  and 
bx  by  2 

y-directions,  respectively,  lb/ft 

S  ,  S  ,  and  S  =  radiation  stresses  which  arise  because  of  the  ex- 
xx  ^  cess  momentum  flux  due  to  waves  (refer  to  Longuet- 

Higgins  and  Stewart  (1964)  for  their  significance) 
lb/ ft 

T  =  lateral  shear  stress  due  to  turbulent  mixing, 

xy  » 

lb/ftz 

Note  that  the  condition  rj  >  0  is  known  as  "setup"  and  r|  <  0  is  called 
"setdown. " 

Bottom  friction 

7.  At  present,  the  numerical  model  uses  a  linear  formulation  for  fric¬ 
tion  (Longuet-Higgins  1970).  Thus 


Tbx  =  PC  "horbP  U 


T,  =  pc  <  u  .  >  V 
by  orb 


where  c  is  a  drag  coefficient  (of  the  order  of  0.01)  and  <|uor{3|>  t^ie 
time  average,  over  one  wave  period,  of  the  absolute  value  of  the  wave  orbital 
velocity  at  the  bottom.  From  linear  wave  theory, 


<  Uorb  >  T  sinh  klhl 


*  A  table  of  factors  for  converting  US  customary  units  of  measurement  to 
metric  (SI)  units  is  presented  on  page  3. 


where 


rvy*  v 


v-..  \~v\  \-  v-  v  v ~-.-_v 


H  =  local  wave  height,  ft 
T  =  wave  period,  sec 
k  =  local  wave  number,  2n/L  ,  1/ft 
|h|  =  local  still-water  depth,  ft 

Equations  4  and  5  are  based  on  the  assumption  that  the  velocity  components  U 
and  V  of  the  current  are  small  compared  with  the  wave  orbital  velocity, 
<|uorb|>  *  In  the  future,  the  numerical  model  can  be  easily  adapted  to  other 
formulations  for  friction  such  as  nonlinear  friction. 

Radiation  stresses 

8.  As  mentioned  previously,  the  radiation  stresses  are  of  major  im¬ 
portance  since  they  furnish  the  main  forces  for  creating  wave-induced  currents. 
Referring  to  Longuet-Higgins  (1970),  for  monochromatic  waves,  they  are  defined 
in  terms  of  the  local  wave  climate  as  follows: 


S  =  E 

?2n  -  I')  cos2  0  +  (n  -  ^  sin2  ©I 

XX 

V  2)  V  2/  J 

S  =  E 

n  cos  0  sin  0 

xy 

S  =  E 

^2n  -  4^  sin2  0  +  (n  -  4 )  cos2  0 

yy 

\  2/  V  2}  J 

E  =  |  PgH2 


and 


(7) 

(8) 

(9) 


(10) 


(ID 


Note  that  n  is  the  ratio  of  wave  group  celerity  to  phase  celerity,  0  is 
the  local  wave  direction  defined  as  shown  in  Figure  1,  and  E  is  the  wave 
energy  density,  lb/ft.  For  the  numerical  model  described  herein,  the  values 
of  H  ,  k  ,  and  0  are  obtained  by  using  a  considerably  modified  form  of  the 
wave  climate  program  developed  by  Noda  et  al.  (1974).  This  particular  program 
has  the  advantage  that  H  ,  k  ,  and  0  can  be  computed  at  the  centers  of  the 
cells  of  a  rectangular  numerical  grid  and  wave  breaking  and  decay  are  accounted 
for  by  a  breaking  index  model  for  wave  heights  in  the  surf  zone. 


..•VVA’.vV.^V. 
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N  S  S 


Lateral  shear 


9.  In  the  numerical  model,  the  coordinate  scheme  is  chosen  such  that  x 
is  positive  in  the  offshore  direction  and  y  is  approximately  in  the  along¬ 
shore  direction.  An  eddy  viscosity  formulation  is  chosen  for  the  lateral 
shear.  The  eddy  viscosity  is  assumed  to  be  anisotropic.  Denoting  and 

£y  as  the  eddy  viscosities  in  x-  and  y-directions ,  respectively,  in  general, 

e  is  assumed  to  be  a  function  of  x  and  y  and  e  a  constant, 

x  y 

Accordingly, 


r  3u,  av 

Cy  9y  £x  dx 


For  plane  beach  applications  with  lateral  mixing,  the  eddy  viscosity  was 

assumed  to  vary  within  the  surf  zone  in  the  manner  suggested  by  Longuet-Higgins 
(1970): 


£x  =  NLHX  ffi" 


where 

=  an  empirical  coefficient  that  varies  in  the  range  0.0  to  0.016 
x  =  distance  from  the  shoreline 

For  locations  offshore  of  the  breaker  line,  was  kept  constant  and  equal 

to  its  value  at  the  breaker  line.  For  field  applications,  the  eddy  viscosity 
ex  was  chosen  according  to  the  following  relation  given  by  Jonsson,  Skovgaard, 
and  Jacobsen  (1974): 


£  =  H2gT 

*  4n2|h| 


COS2  0 


This  represents  twice  the  value  used  by  Thornton  (1970).  It  was  believed  that 
for  field  situations,  Equation  14  represented  the  eddy  viscosities  more  realis¬ 
tically  than  the  relation  (Equation  13)  for  plane  beaches.  The  value  of 
was,  in  general,  taken  to  be  equal  to  the  value  of  at  the  deepest  part 

(usually  near  the  offshore  boundary)  of  the  numerical  grid.  The  numerical 
model  is  flexible  enough  to  permit  other  formulations  for  eddy  viscosity  in 
the  future,  as  our  understanding  improves. 

Variable  grid 

10.  One  major  advantage  of  WIFM  and  therefore  CURRENT  is  that  the  size 
of  the  grid  cells  may  be  varied  smoothly  in  both  horizontal  directions.  Thus 
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the  grid  may  be  made  finer  in  regions  of  greater  interest  such  as  the  surf 
zone,  inlets,  etc.,  and  coarser  in  regions  of  less  importance.  For  this  pur¬ 
pose,  a  piecewise  reversible  transformation  is  used.  The  mapping  for  each 
coordinate  direction  is  independent  of  the  other.  For  mapping  from  the  real 
or  physical  space  (x,  y)  to  the  computational  space  (a^,  a.,),  a  mapping  func¬ 
tion  of  the  form 


x  =  at  +  bjffj  (15) 

is  used.  A  similar  mapping  function  is  employed  for  the  y-direction.  Here 
a^  ,  bj  ,  and  c^  are  coefficients  whose  values  change  from  region  to  region. 
The  mapping  transforms  the  variable  grid  in  real  space  to  a  uniform  grid  in 
computational  space.  The  transformation  is  such  that  all  derivatives  are 
centered  in  a-space.  During  the  mapping  process,  both  the  variable  and  its 
derivative  must  match  at  the  common  boundary  between  two  regions.  The  mapping 
is  actually  accomplished  by  an  iterative  procedure,  using  an  interactive  pro¬ 
gram  called  MAPIT  developed  at  WES.  Figure  2  shows  an  example  of  the  variable 
grid. 


REAL  SPACE 


1  A&2 

COMPUTATIONAL  SPACE 


MAPPING  FUNCTIONS 

x  *  «1  +  b,  of1 

y  -  a2  +  b2  «22 

Figure  2.  An  example  of  variable  grid 

11.  By  using  the  mapping  function  in  the  x-  and  y-directions ,  the  equa¬ 
tions  of  motion  in  computational  space  may  be  written  as: 


i  **•  *"■  **•  •'■**«*  v’  •/  «,*  •'  O  ».*  », 


/ 


Momentum 


V  +  —  UV  +  —  W  +  &-  rj  +  K\—  (s  )  +  —  (s  \  1 

t  Ml  orj  m2  «2  m2  a2  Pd  by  pd  Mj  V  M2Vyy/a| 

-e  —  —  V  +  ( — \  V  -£  —  —  U  -  ^  /e  N  v  =  0 

x  »*i  »*i  Vl  \V  al  y  M1  ^2  V2  M|  '  '<*.  "l 
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Continuity 


n  +  —  +  —  (vA\  =  0 

1  Mj  V  hl  m2  V  )a2 


where  the  subscripts  t  ,  ,  and  or 2  indicate  partial  derivatives  with 

respect  to  time,  and  a 2  ,  respectively,  and 

dx  .  Cl_1  . 

Mi  =  §5;  =  bici“i  ( 


9y  c2-1 

M2  =  =  b2C2°2 


Variables  and  p2  are  expansion  coefficients  connected  with  the  stretch¬ 

ing  of  the  uniform  computational  grid  to  the  variable  grid  in  real  space. 

Note  that  in  obtaining  Equations  16  and  17,  the  assumptions  made  in  paragraph  9 
were  used. 

12.  The  nonlinear  advective  terms  in  the  equations  of  motion  often  pose 
stability  problems.  These  terms  are  handled  in  the  present  model  by  using  the 
Stabilizing  Correction  (SC)  scheme.  This  scheme  will  be  described  in  a  later 


section.  The  eddy  viscosity  terms  also  can  cause  difficulties  during  the 
numerical  computation.  The  finite  difference  scheme  selected  and  the  formula¬ 
tion  for  eddy  viscosity  adopted  in  the  model  will  minimize  such  difficulties 
and  stability  problems.  However,  the  user  must  exercise  caution  and  judgment 
in  selecting  appropriate  time-  and  space-steps,  and  eddy  viscosity  coefficients 
for  the  phenomena  being  simulated. 


PART  III:  COMPUTATIONAL  TECHNIQUES 

Implicit  Method 

13.  In  order  to  solve  the  problem  under  consideration  on  a  digital  com¬ 
puter,  the  differential  equations  (Equations  16-18)  have  to  be  expressed  in  a 
finite  difference  form.  In  the  present  case,  an  alternating  direction,  im¬ 
plicit,  finite  difference  scheme  is  employed.  In  view  of  the  presence  of  the 
nonlinear  advective  terms,  a  particular  implicit  scheme  known  as  the  SC  scheme 
is  used.  The  basic  idea  of  the  scheme  is  as  follows.  The  time  level  is  indi¬ 
cated  by  a  superscript  k  .  The  scheme  involves  variables  at  three  time 
levels.  The  values  of  the  variable  at  time  levels  k-1  and  k  are  known 
from  previous  computations  or  prescribed  initial  conditions.  To  advance  the 
solution  from  time  level  k  to  the  new  time  level  k+1  ,  we  introduce  an 
intermediate  time-level  solution  denoted  by  the  superscript  *  .  We  operate 
on  Equations  16-18  in  a  two-step  procedure.  In  the  first  step,  we  sweep  the 
rectangular  grid  in  the  x(a^)  direction,  advancing  the  solution  from  time 
level  k  to  *  .  Next,  we  sweep  the  grid  in  the  yCo^)  direction,  advancing 
the  solution  from  time  level  *  to  k+1  .  The  two  sweeps  together  constitute 
a  full  time-step,  At  . 


Double-Sweep  Technique 


14.  Before  going  into  the  deta 
notation  used  for  individual  cells  of 
Let  Ax  and  Ay  denote  the  cell  di¬ 
mensions  in  real  space  in  the  x-  and 
y-directions ,  respectively.  These 
dimensions  may  vary  from  cell  to  cell 
Let  the  corresponding  dimensions  in 
computational  space  be  Aa^  and 
A«2  .  These  dimensions  are  the  same 
for  all  the  cells  in  the  grid.  Let 
m  and  n  denote  indices  correspond¬ 
ing  to  the  center  of  an  arbitrary 
cell  (Figure  3).  All  the  variables 


Is  of  the  double  sweep  technique,  the 
the  rectangular  grid  will  be  defined. 


Figure  3.  Cell  notation 


except  the  velocities  U  and  V  are  defined  at  the  cell  centers.  Velocities 

U  and  V  are  defined  at  cell  faces  m+1/2  and  n+1/2  ,  respectively.  In 

the  x-sweep,  the  x-momentum  equation  is  centered  about  the  cell  face  m+1/2 

and  the  continuity  equation  is  centered  about  the  center  of  the  cell  and  the 

k+1 

two  equations  are  solved,  using  in  the  process  the  result  U*  =  U  .At  the 

_  fc+i 

end  of  this  sweep,  q*  and  U  are  known.  Next,  we  sweep  the  grid  in  the 

y-direction.  In  this  sweep,  the  y-momentum  equation  is  centered  about  the  cell 

face  n+1/2  and  the  continuity  equation  about  the  cell  center.  Upon  solving 

-k+1  k+1 

the  two  equations,  q  and  V  for  each  cell  are  obtained.  Thus  the  two 
sweeps  together  complete  the  solution  for  qk+*  ,  Uk+1  ,  and  Vk+1  . 

15.  Even  though  the  SC  scheme  has  been  described  so  far  in  terms  of  the 
(x,  y)  coordinate  system  for  convenience,  in  reality  we  have  to  apply  the  tech¬ 
nique  to  the  equations  of  motion  in  the  computational  (a^,  a^)  space.  After 
the  application  of  the  technique,  the  following  finite  difference  equations 
result.  For  the  rest  of  PART  III,  we  will  drop  the  bar  over  q  for  conve¬ 
nience.  For  the  a^(x)-sweep  (taken  along  a  grid  cell  column  parallel  to  the 
a^-axis),  we  have 

1  ,„k+l  TTk-l  Uk  ,  ,„k.  ?k  *  n,kx 

2At  (U  u  )  +  2^4^  62aj(U  J  2M2Aa2  62a2(U  ) 
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A„  =  r1k'1  -  iT£T  6*  (uk_1  ^  -  fTST  6C  (vk_1  *k)  at  (">■“>  <30> 

H1  1  h22  2 

Vk  =  —  /yk  +  Vk  +  Vk  +  Vk  \  f31) 

n, m+1/2  4  ^  n-l/2,m  n+l/2,m  n-l/2,m+l  n+l/2,m+l) 

17.  Consider  the  set  of  cells  for  which  the  index  n  is  constant  and 
equal  to  N  .  Suppose  at  the  upper  boundary  cell  (m  =  M) ,  the  velocity 
M+1/2  *s  always  known.  Similarly  suppose  at  the  lower  boundary  cell 
(m  =  L)  the  water  level  rtT  T  is  always  known.  Then  the  set  of  equations  for 

N  ,  L 

all  the  cells  may  be  written  in  the  following  matrix  form  if  we  drop  the  com¬ 
mon  subscript  N  : 


aM+l/2  aM+l 
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18.  Since  the  first  matrix  on  the  left-hand  side  of  Equation  32  is  tri¬ 
diagonal,  the  above  matrix  equation  can  be  solved  by  recursion.  In  general, 
the  recursion  relations  may  be  written  as 


m  m+1/2 
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19.  Since  in  the  FORTRAN  computer  language,  fractional  indices  are  not 
possible,  a  new  integer  index  system  is  adopted  in  the  program.  Thus  all  the 
variables  defined  at  the  center  and  the  faces  m+1/2  and  n+1/2  of  a  cell 
(n,m)  will  be  designated  by  the  integer  indices  (N,M).  The  only  exceptions 
are  the  expansion  coefficients  p^  and  M2  which  are  defined  at  cell  centers 
and  faces.  For  these,  the  following  index  system  is  adopted.  For  example, 

Mj  at  the  center  of  cell  (n,m)  is  designated  by  the  index  2M-1  ,  whereas  Pj 
at  the  face  m+1/2  is  denoted  by  the  index  2M  ,  and  similarly  for  P2  . 

Using  this  new  notation,  the  expanded  form  of  the  recursion  coefficients  for 
the  a  (x)-sweep  may  be  written  as  follows: 


Using  the  same  notation,  the  solution  (Equations  33  and  34)  may  be  written  as 


nN,M  “PMUN,M  +  QM 


K  *  1  ^ 

UN,M-1  =  _RM-inN,M  +  SM-1 


For  any  given  N  ,  the  recursion  coefficients  P  ,  Q  ,  R  ,  and  S  are  com¬ 
puted,  using  Equations  36-41,  in  succession  between  the  boundaries  in  the  di¬ 

rection  of  increasing  or^(x)  •  The  values  of  these  coefficients  at  the  bound¬ 
aries  depend  on  the  types  of  boundary  conditions  encountered.  Once  all  the 

k+1 

coefficients  for  a  given  N  have  been  determined,  the  values  of  q*  and  U 
for  all  the  cells  in  the  column  are  computed,  using  Equations  41  and  42,  in 
the  direction  of  decreasing  a^(x)  •  We  next  go  to  the  next  higher  value  of 

N  ,  and  so  on  until  the  whole  grid  is  swept  in  the  a^(x)-direction. 

20.  The  development  of  the  finite  difference  equations  and  the  recursion 
relations  for  the  a^yj-sweep  is  similar  to  that  for  the  a^(x)-sweep.  In 
this  case,  using  the  same  notation  as  before,  the  recursion  coefficients  may 
be  written  as 
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The  corresponding  solution  may  be  expressed  as 


nk+1  =  -P  vk+1  +  Q 
NPM  111,(1 


Initial  and  Boundary  Conditions 


21.  In  order  to  solve  the  problem  under  consideration,  appropriate 
initial  and  boundary  conditions  must  be  specified.  For  the  examples  reported 
here,  an  initial  condition  of  rest  was  chosen  so  that  r|  ,  U  ,  and  V  are 
zero  at  the  start  of  the  calculations.  To  avoid  shock,  the  radiation  stress 
gradients  were  gradually  built  up  to  their  full  values  over  a  number  of  time- 
steps.  The  numerical  computation  was  stopped  when  a  steady  state  was  deemed 
to  have  been  reached. 


22.  The  numerical  model  permits  different  types  of  boundary  conditions; 
among  these  are  the  following: 

a.  "No  flow"  (wall).  This  type  of  boundary  condition  is  used  at 
closed  boundaries  such  as  the  still-water  line  on  beach  and  at 
impermeable  structures.  The  normal  velocity  is  set  to  zero  in 
this  case. 

b.  Discharge .  This  type  of  condition  is  used  by  WIFM  at  open 
boundaries.  The  variation  of  discharge  with  time  is  prescribed 
along  boundary  cells.  For  the  wave-induced  current  model,  this 
condition  is  never  used  since  the  discharge  is  an  unknown  even 
at  the  boundaries. 

c.  Elevation  (tide).  This  type  of  condition  is  used  by  WIFM  at 
open  boundaries.  The  variation  of  the  surface  elevation  with 
time  is  prescribed  along  boundary  cells.  For  the  wave-induced 
current  model,  this  condition  is  never  used  except  possibly  at 
the  offshore  boundary  since  the  setup  is  an  unknown  quantity. 
Even  at  the  offshore  boundary,  a  radiation  boundary  condition 
was  found  to  be  preferable. 

d.  Uniform  flux.  In  this  type  of  open  boundary  condition,  the 
flux  at  a  boundary  cell  is  made  equal  to  that  at  the  next  in¬ 
terior  cell.  Thus  the  condition  assumes  3(Ud)/3x  =0  or 
3(Vd)/3y  =0  at  the  boundary.  This  type  of  condition  is  used 
for  the  lateral  boundaries  since  it  is  a  passive  condition. 

e.  Radiation.  This  open  boundary  condition  requires  that  any 
transients  developed  initially  inside  the  numerical  grid  should 
propagate  out  of  the  grid  as  gravity  waves.  It  is  of  the  form 
3q/3t  +  C(3r)/3x)  =  0  where  C  is  the  phase  speed  of  a  surface 
disturbance  r|(x,t)  •  It  is  often  used  by  the  wave-induced  cur¬ 
rent  at  the  offshore  boundary  and  is  found  preferable  to  a  wall 
or  constant  elevation  condition  there.  Both  of  the  latter  con¬ 
ditions  are  highly  reflective  and  as  a  result  the  transients 
tend  to  bounce  back  and  forth  between  the  offshore  and  nearshore 
boundaries  and  take  a  long  time  to  damp  out.  On  the  other  hand, 
the  radiation  condition  seems  to  work  quite  well  and  allows  the 
transients  to  propagate  out  of  the  grid  which  permits  the  set- 
down  at  the  offshore  boundary  to  assume  an  appropriate  value. 

23.  The  boundary  conditions  frequently  used  in  the  wave-induced  current 
model  are  illustrated  in  Figure  4. 

24.  At  present,  the  model  allows  for  subgrid  barriers  such  as  jetties, 
provided  they  are  impermeable  and  nonovertopping.  The  program  essentially 
sets  to  zero  the  velocity  component  normal  to  the  appropriate  cell  face. 


Figure  4.  Boundary  conditions  used  in  numerical 
model  CURRENT 
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PART  IV:  VALIDATION  OF  MODEL 

Tests  for  Idealized  Conditions 


25.  To  develop  confidence  in  the  validity  of  the  model  and  the  accuracy 
of  its  results,  several  tests  were  run  on  the  model  and  comparisons  were  made 
between  model  results  and  available  laboratory  data  and  analytic  solutions. 

All  of  these  tests  were  for  plane  beaches,  for  which  the  coordinate  scheme  is 
chosen  such  that  the  y-axis  coincides  with  the  still-water  line  in  beach  and 
the  x-coordinate  is  measured  from  the  still-water  line  (Figure  5).  Note  that 
for  plane  beaches,  there  is  no  variation  in  the  alongshore  (y)  direction. 

Plane  beach:  normal  incidence 

26.  The  model  was  run  for  a  case  of  normal  incidence  on  a  plane  smooth 
laboratory  beach,  reported  by  Bowen,  Inman,  and  Simmons  (1968).  The  conditions 
were  as  follows:  T  =  1.14  sec,  deepwater  wave  height  Hq  =  6.45  cm,  and  beach 
slope  s  =  1:12  .  To  run  this  case  on  the  model,  a  variable  rectangular  grid 
with  overall  dimensions  of  approximately  40  m  by  30  cm  (the  laboratory  channel 
was  40  m  long)  was  used  with  Ao^  =  Aof^  =  10  cm  and  At  =  0.05  sec.  The  grid 
was  3  cells  wide  in  the  alongshore  direction  and  50  cells  long  in  the  offshore 
direction.  In  this  example,  walls  were  used  for  the  lateral  boundaries  as  well 
as  the  offshore  boundary  to  correspond  to  the  laboratory  situation.  Since  for 
normal  incidence,  the  velocities  U  and  V  would  be  zero  everywhere  corre¬ 
sponding  to  the  steady  state,  advection,  eddy  viscosity,  and  friction  terms  were 
turned  off  in  the  model.  The  solution  allowed  for  the  effect  of  setup  on  the 
wave  heights  in  the  surf  zone.  As  the  solution  proceeded,  since  n  changed, 
the  wave  heights  for  cells  in  the  surf  zone  were  computed  afresh  for  each  time- 
step  by  using  H  =  y(  |  h|  +  r|),  where  y  is  a  breaking  index  and  the  radiation 
stresses  were  changed  accordingly.  As  suggested  by  Bowen,  Inman,  and  Simmons 
(1968),  a  y  of  1.15  was  used.  A  buildup  time  of  10  At  was  used  at  the 
start.  A  comparison  of  the  steady-state  setup  values  from  the  model  (after 

150  At)  with  those  observed  by  Bowen,  Inman,  and  Simmons  (1968)  is  shown  in 
Figure  6.  As  shown,  there  is  excellent  agreement  in  the  offshore  region.  In 
the  surf  zone,  the  numerical  model  predicts  higher  setups  than  those  observed. 
This  is  not  surprising  since  the  model  does  not  allow  flooding  and  runup.  It 
is  to  be  noted  that  the  slope  of  the  mean  waterline  in  the  surf  zone  is  ap¬ 
proximately  the  same  in  both  cases. 
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Figure  5.  Definition  sketch  for  a  plane  beach:  cross  section  and  plan 
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Figure  6.  Comparison  of  numerical  solution  for 
setup  with  experimental  data 


Plane  beach:  oblique  incidence 

27.  For  this  case,  a  plane  beach  of  constant  bottom  slope  s  =  1:30 
was  selected.  A  monochromatic  wave  with  the  following  deepwater  characteris¬ 
tics  was  chosen:  T  =  12  sec,  Hq  =  10  ft,  and  angle  of  incidence  in  deep 
water,  6^  =  20  deg.  A  drag  coefficient  c  of  0.01  and  a  breaking  index  y 
of  0.82  were  used  in  the  model.  A  uniform  grid  with  Ax  =  Ay  =  60  ft  was 
used  for  most  of  the  runs.  It  was  6  cells  wide  in  the  alongshore  direction 
and  100  cells  long  in  the  offshore  direction.  Uniform  flux  and  radiation 
boundary  conditions  were  used  for  the  lateral  and  offshore  boundaries,  re¬ 
spectively.  The  buildup  time  varied  from  15  to  50  At  ,  depending  on  the 
time-step  At  used. 

28.  First  the  model  was  run  without  allowing  for  the  effect  of  setup  on 
wave  heights  and  radiation  stresses.  Mixing  and  advection  were  ignored.  A 
time-step  of  0.5  sec  was  used.  The  steady-state  velocity  distribution  obtained 
(after  800  At)  is  compared  with  the  triangular  distribution  of  Longuet-Higgins 
(1970)  in  Figure  7.  There  is  good  agreement.  Note  that  for  positive  0  ,  V 
will  be  negative  for  our  coordinate  scheme.  Later  a  finer  grid  (Ax  =  Ay 

=  30  ft)  with  a  At  of  0.25  sec  was  used.  As  shown  in  Figure  7,  as  the  grid 
is  made  finer,  the  numerical  solution  tends  to  approach  the  analytic  solution. 

29.  The  effect  of  setup  was  taken  into  account  next.  A  time-step  of 


26 


Figure  7.  Plane  beach:  solution  for  longshore  current 
neglecting  the  effect  of  setup 

1.5  sec  was  used  for  this  case.  The  velocity  distribution  from  the  model  is 
compared  with  the  corresponding  analytic  solution  in  Figure  8.  There  is  good 
agreement.  Note  that  the  numerical  solution  goes  to  zero  at  the  still-water 
line  because  a  wall  was  assumed  there.  On  the  other  hand,  the  Longuet-Higgins 
(1970)  solution  goes  to  zero  at  the  setup  line.  To  plot  his  solution,  the  dis¬ 
tance  from  the  still-water  line  to  the  setup  line  was  estimated  by  using  a  re¬ 
lation  provided  by  Dalrymple,  Eubanks,  and  Birkemeier  (1977). 

30.  The  effect  of  lateral  mixing  was  studied  next,  without  taking  the 
effect  of  setup  into  account.  A  time-step  of  5.0  sec  was  used  for  these  runs. 
The  mixing  parameter  P  of  Longuet-Higgins  (1970)  was  varied  between  0.01  and 
0.4.  Note  that  P  is  defined  as 


Y  c 


(52) 
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-  ANALYTICAL  SOLUTION  (LONGUET- HIGGINS  1970) 

- COARSE  GRID  Ax  -  Ay  =  60  FT,  At  =  1.5  SEC 
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Figure  8.  Plane  beach:  solution  for  longshore  current 
considering  the  effect  of  setup 

Figure  9  shows  the  effect  of  P  on  the  numerical  solution.  As  expected,  the 
magnitude  of  the  peak  decreases,  the  peak  moves  closer  to  the  shoreline,  and 
the  velocities  offshore  of  the  breaker  line  increase  as  P  increases. 


Difficulties  Involved  in  Application  to  Field  Situations 

31.  While  it  is  relatively  easy  to  apply  a  numerical  current  model  to 
idealized  cases,  one  must  face  several  difficulties  in  applying  the  model  to 
field  situations.  Among  these  is  the  highly  irregular  nature  of  the  bathym¬ 
etry,  especially  near  inlets  where  channels  and  shoals  exist.  The  topography 
must  be  smoothed  to  a  certain  extent  in  order  for  the  wave  climate  and  wave- 
induced  current  models  to  work  properly;  yet,  one  must  be  careful  not  to  com¬ 
pletely  change  the  basic  features  of  the  topography.  The  shoreline  as  well  as 
the  breaker  line  may  be  irregular  and  may  be  oblique  to  the  grid  axes.  There 
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Figure  9.  Plane  beach:  effect  of  mixing  parameter  P  on  the 

numerical  solution 

may  be  more  than  one  breaker  line.  There  are  problems  connected  with  discreti 
zation  of  the  shoreline  and  breaker  line(s).  Selection  of  appropriate  values 
for  empirical  coefficients  such  as  friction  and  eddy  viscosity  coefficients 
and  breaking  index  is  not  easy.  There  are  problems  in  connection  with  the 
wave  climate  model  also,  especially  if  wave-current  interactions  are  to  be 
taken  into  account. 

A  Particular  Field  Application 

32.  In  order  to  demonstrate  the  applicability  of  the  numerical  model  to 
field  situations,  the  case  of  Oregon  Inlet,  North  Carolina,  was  selected. 
Oregon  Inlet  is  a  tidal  inlet  in  a  barrier  island  system.  Behind  the  inlet 
toward  the  mainland  is  Pamlico  Sound.  Most  of  the  problems  mentioned  in  the 
previous  paragraph  had  to  be  addressed  and  solved  satisfactorily  in  this 


application.  For  purposes  of  the  numerical  simulation,  a  rectangular  region 
approximately  62,400  ft  long  in  the  alongshore  direction  and  29,400  ft  wide  in 
the  offshore  direction  was  considered.  It  included  a  portion  of  Pamlico  Sound. 
The  variable  grid  used  for  the  simulation  is  shown  in  Figure  10.  The  grid  was 
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Figure  10.  Numerical  grid  used  for  Oregon  Inlet, 
North  Carolina,  simulation 


77  cells  long  in  the  alongshore  direction  and  54  cells  wide  in  the  offshore 
direction.  It  may  be  noted  that  the  minimum  cell  widths  in  the  alongshore  and 
offshore  directions  were  400  and  100  ft,  respectively.  These  widths  were  used 
near  the  inlet  and  surf  zone,  respectively.  Note  that  Adfj  =  At^  =  100  ft. 

The  topography  used  in  the  simulation  corresponding  to  this  grid  is  shown  in 
Figure  11.  The  elevations  are  shown  in  feet  and  the  datum  is  mean  low  water 
(mlw) .  Several  points  must  be  mentioned  about  this  three-dimensional  per¬ 
spective  plot.  First,  the  vertical  dimensions  are  highly  exaggerated  compared 
with  the  horizontal;  secondly,  the  elevations  are  plotted  in  the  computational 
space  and  not  the  physical  space--so  the  horizontal  dimensions  are  distorted. 
The  topography  was  somewhat  modified  compared  with  the  actual  topography,  with 
respect  to  the  depths  near  the  offshore  boundary  and  the  land  elevations  on 
the  islands.  In  spite  of  these  factors,  Figure  11  helps  one  to  visualize  the 
irregular  nature  of  the  bathymetry.  Also,  the  locations  of  the  channels  and 
shoals  in  the  region  of  the  inlet  are  shown  clearly  in  the  figure. 

33.  A  monochromatic  wave  with  a  height  of  11.39  ft,  period  of  8.0  sec, 
and  6  =  51.1  deg  in  60-ft  depth  of  water  was  selected  for  the  simulation 
(the  depth  of  water  at  the  offshore  boundary  of  the  numerical  grid  was  60  ft). 
This  wave  corresponded  to  the  significant  wave  during  a  part  of  the  Ash  Wednes¬ 
day  storm  of  March  1962  at  the  inlet.  In  this  case,  besides  using  "no  flow" 
conditions  at  the  shoreline,  a  radiation  boundary  condition  offshore,  and  uni¬ 
form  flux  boundary  conditions  at  the  lateral  boundaries,  a  uniform  flux  bound¬ 
ary  condition  was  used  over  a  part  of  the  inland  side  of  the  sound,  while  the 
rest  of  the  sound  was  closed  off.  A  time-step  At  of  18.0  sec  and  a  drag 
coefficient  c  of  0.01  were  used  in  the  numerical  model.  The  breaking  index 
Y  was  chosen  according  to  the  breaking  criterion  employed  by  Noda  (1974): 


where  L  corresponds  to  the  wavelength  and  the  subscript  b  indicates  values 
at  breaking.  A  buildup  time  of  15  At  was  used  at  the  start.  The  eddy  vis¬ 
cosity  e  was  chosen  according  to  Equation  14  and  the  eddy  viscosity  £  was 
x  y 

set  equal  to  the  value  of  £x  at  the  offshore  boundary.  For  the  case  under 
consideration,  the  complete  equations  (Equations  1,  2,  and  3)  were  solved.  An 
approximate  steady  state  was  reached  after  67  At  .  Figures  12  and  13  represent 


VELOCITY  =  5  FT/SEC 


Figure  13.  Velocity  vector  plot  for  Oregon  Inlet, 

North  Carolina,  simulation 

the  corresponding  mean  water  levels  and  velocity  vectors  plotted  on  the  grid 
in  the  computational  space.  The  velocity  vectors  are  plotted  for  every  other 
cell  in  each  coordinate  direction.  To  avoid  confusion,  the  plotting  of 
velocities  with  magnitudes  less  than  0.1  ft/sec  is  suppressed. 

34.  Referring  to  Figures  11,  12,  and  13,  let  us  first  consider  the  two 
portions  of  the  beach  away  from  the  inlet.  The  shorelines  in  these  regions 
are  approximately  straight  and  the  contours  are  approximately  straight  and 
parallel.  As  we  approach  the  shoreline  from  offshore,  there  is  a  small  setdown 
followed  by  a  setup.  The  velocities  are  mainly  alongshore  and  the  velocity 
distribution  is  similar  to  that  for  a  plane  beach  except  that  it  exhibits  two 
peaks  at  some  locations. 

35.  The  situation  is  more  complicated  in  the  region  of  the  inlet  (the 
central  part  of  the  grid).  Here  the  breaker  line  is  farther  offshore.  The 


depth  in  the  main  channel  decreases  first  and  increases  later  as  we  go  toward 
the  inlet.  Because  of  these  factors,  the  water  sets  up  around  the  inlet  and 
tends  to  create  a  flow  into  the  inlet  through  the  various  channels,  as  one 
would  naturally  expect.  A  part  of  the  main  alongshore  flow  goes  around  the 
channels  and  shoals  to  the  other  side. 

36.  Near  the  shoals,  the  patterns  of  mean  water  level  and  velocity  are 
irregular.  This  is  because  the  waves  refract  around  the  shoals  and  break, 
creating  locally  setups  and  currents  that  do  not  necessarily  conform  to  the 
general  patterns.  As  the  waves  go  toward  the  islands,  they  re-form  because 
the  depth  increases. 

37.  Figures  12  and  13  do  not  reflect  the  influence  of  tides  and  fresh¬ 
water  flows  through  the  inlet.  In  nature,  these  phenomena  tend  to  modify  the 
patterns  shown  in  these  figures. 


Summary 


38.  The  various  tests  for  idealized  situations  and  the  successful  com¬ 
parisons  to  analytic  solutions  and  experimental  data  indicate  that  the  numeri¬ 
cal  model  "behaves  properly"  and  yields  valid  and  accurate  results.  In  the 
case  of  the  field  situation  for  Oregon  Inlet,  there  are  no  field  data  avail¬ 
able  with  which  the  model  results  can  be  compared.  However,  for  this  very 
complicated  case,  the  model  yields  results  that  appear  to  be  reasonable. 


PART  V:  MODEL  INPUT 


General  Description 

39.  The  input  data  for  CURRENT  are  discussed  in  this  part.  The  model 
has  been  tested  using  the  foot-pound-second  system  of  units  only.  The  number 
of  time-steps  computed  is  indicated  by  an  integer  variable  ITIME.  The  grid 
index  system  is  set  up  with  N  approximately  in  the  alongshore  direction  and 
M  increasing  in  the  offshore  direction.  If  a  variable  grid  is  used,  the  ex¬ 
pansion  coefficients  must  be  defined  for  both  cell  centers  and  cell  faces  and 
are  obtained  from  the  programs  MAPIT  and  GRID  developed  at  WES. 

Setting  Matrix  Dimensions 

40.  In  the  CRAY  computer  system  on  which  the  model  was  developed,  matrix 
dimensions  are  changed  for  a  particular  application  by  using  PARAMETER  cards. 
The  following  parameters  must  always  be  set.  This  is  done  by  using  the  UPDATE 
feature  of  CRAY.  Using  text  editor,  the  following  statements  in  "ID  PARMRUN 

in  update  file  UPDL2RV  should  be  changed  appropriately  for  the  application: 

PARAMETER (NDIM=  ,  MDIM=  ,  LDIM=  ,  IFRC=  ) 

PARAMETER ( I PUT=  ,  ITDS=  ,  JPUT=  ,  JFLS=  ) 

PARAMETER (NGAGES=  ,  NBCELS=  ,  NBARCL=  ) 

The  significance  of  the  various  parameters  is  explained  below.  Normally,  only 
the  parameters  marked  with  an  asterisk  have  to  be  changed.  The  rest  may  be 
left  at  the  values  given  with  the  program. 

NDIM*  =  Grid  dimension  (number  of  cells)  in  y-direction 
MDIM*  =  Grid  dimension  (number  of  cells)  in  x-direction 
LDIM*  =  Larger  of  NDIM  and  MDIM 

IFRC*  =  Dimension  of  friction  array  for  Chezy  values  (>4*DMAXG, 
where  DMAXG  is  the  maximum  total  water  depth  (ft)  that 
will  be  experienced  anywhere  in  the  grid  during  the 
numerical  simulation) 

IPUT*  =  Dimension  of  elevation  (tidal)  forcing  array-array  must 
include  a  value  for  each  time-step 

ITDS*  =  Number  of  different  elevation  (tidal)  forcing  arrays 
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JPUT 

JFLS 

NGAGES 

NBCELS* 

NBARCL* 


Dimension  of  discharge  forcing  arrays 

Number  of  different  discharge  forcing  arrays 

Number  of  special  gage  locations  for  which  results  are 
to  be  printed 

Number  of  boundary  cells  (where  boundary  conditions  are 
to  be  applied) 

Number  of  barrier  cell  faces 


Input  Data 


41.  The  input  data  for  CURRENT  have  been  assembled  into  card  groups 
that  may  consist  of  one  or  more  data  cards.  Some  groups  are  optional  and  thus 
each  group  is  marked  with  an  R  (required)  or  0  (optional).  The  following  text 
indicates  for  each  input  card  group  the  necessity  code,  FORTRAN  format,  the 
variables  involved,  and  a  brief  description  of  the  variables. 


Necessity  Card  Group 

Code _  (Format)  Variable 


Description 


R  1A 

(15) 


NDTAP  Number  of  file  from  which  input  data  are 

to  be  read  in--usually  set  to  95 


NOTE:  All  the  following  card  groups  are  on  a  separate  file  defined 
by  NDTAP  (if  NDTAP^5).  Group  1A  must  be  on  the  system  input 
file  (File  5,  i.e.  card  reader) 


R 

R 


R 


IB 

(8A8) 

1C 

(NAMELIST 

$WAVE) 


ITL 

PER 

HTO 

THETAO 


DEPMAX 


GAMMA 


2  NMAX 

(315) 

MMAX 


Identification  title  for  run 


Period  of  the  wave 


Wave  height  at  the  offshore  boundary- 
used  for  calculating  e 

y 

Angle  of  incidence  of  waves  (in  degrees) 
at  the  offshore  boundary,  measured  with 
respect  to  the  x-axis--used  for  calcu¬ 
lating  £y 


Depth  at  the  offshore  boundary  (ft)-- 
used  for  calculating  £ 


Breaking  index  (ratio  of  wave  height  to 
depth  in  surf  zone) 


Grid  dimension  (number  of  cells)  in 
y-direction 

Grid  dimension  (number  of  cells)  in 
x-direction 
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Necessity  Card  Group 

Code  (Format)  Variable 


Description 


R 


R 


3 

ITID 

Number  of  entries  in  input  elevation 

(1615) 

(tidal)  table 

JTID 

Number  of  time-steps  between  entries  in 
input  elevation  (tidal)  table 

NTID 

Number  of  distinct  elevation  (tidal) 
forcing  functions 

MPR 

Print  control  for  initial  conditions 

1- -print  flag  arrays  only 

2 —  print  depth  also 

<0--print  in  addition  flood,  barrier,  and 
elevation  (tidal)  data 

MSURF 

Print  elevation  (tidal)  forcing  function 
in  steps  of  MSURF 

U 

TAU 

Time-step,  At  (sec)  for  one  full  cycle 

(8E10.1) 

DX 

Cell  dimension  Aoi^  (ft) 

DY 

Cell  dimension  Aa  (ft) 

G 

^  2 

Acceleration  due  to  gravity  (ft/sec  ) 

EPSD 

Minimum  amount  of  water  defining  a  dry 
cell  (ft) 

DCON1 

Add  DC0N1  to  water  cell  bed  elevations  to 
translate  datum  (ft) 

DMPX 

Value  of  still-water  land  elevation  as¬ 
signed  artificially  to  areas  that  will 
never  flood  (ft) 

VIS 

A  multiplier  for  eddy  viscosity  terms 
(can  be  set  between  0  and  1) 

XLAND* 

A  value  of  bed  elevation,  h  >XLAND  de¬ 
fines  a  cell  that  will  never  flood  (ft) 
(positive) 

XSCOUR* 

A  value  of  h  <  XSCOUR  defines  a  cell 
that  will  never  go  dry  (ft)  (negative) 

SMAX 

If  n  >  SMAX  ,  stop  computations  and 
print  r|  ;  caution:  set  SMAX  higher  than 
highest  land  cell  elevation 

DMAXG 

Positive  bound  on  maximum  total  water 

depth  that  will  be  experienced  during 
simulation  (ft) — used  for  setting  up  size 
of  friction  matrix 


*  Caution:  Select  XLAND<  bed  elevation  of  lowest  land  cell  and  XSCOUR 

<  depth  |h|  of  shallowest  water  cell  if  flooding/drying  is  not  desirable. 
The  flood/dry  capabilities  of  the  program  have  not  been  tested. 
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Necessity 

Code 


R 

(Continued) 


R 


R 


0 


R 


Card  Group 

(Format)  Variable  _ Description 


4  DC0N2 

(8E10.1) 

(Continued)  DLImIT 


Add  DC0N2  to  elevation  (tidal)  input 
values  to  correspond  to  model  datum  (ft) 

Artificial  cutoff  value  on  bed  elevation 
(h)  for  water  cells  (negative)  (ft) 


5 

(1615) 


MAXTIM 

INTAP 

IDELAY 

IXPAN 

NGAGE 

NFREQ 

NZP 


Number  of  time-steps  to  run  simulation 

=m--save  r|  ,  U  and  V  on  file  1 
every  mAt 

=- 1 --no  data  are  saved  on  file  1 

Delay  saving  information  on  file  1  until 
ITIME=  IDELAY  (note  ITIME  counts  the  num¬ 
ber  of  cycles) 

=0  uniform  grid  in  real  space 
^0  variable  grid--read  in  expansion 
coefficients 

Number  of  special  gage  points  (up  to  20 
gages  are  permitted) 

Information  will  be  printed  at  gage  points 
every  NFREQ  At's--never  set  equal  to  zero! 

Number  of  corrections  to  input  depth 
matrix 


7 

(1615) 


NPRINT  Time-step  index  to  print  hydrodynamics 

(0  *  U  ,  and  V)  over  the  whole  grid — up 
to  32  printouts  are  allowed--two  cards 
must  be  included 


8 

(1615) 


NPOT  y-indices  for  special  gage  points  (NGAGE 

in  number) 

MPOT  x-indices  for  special  gage  points  (NGAGE 

in  number) — start  on  a  new  card 


Note:  Group  8  is  omitted  if  NGAGE=0 


Mannings 's  n  for  each  code  i,  i=l,20-- 
used  for  defining  friction.  Note:  code  1 
is  used  for  all  water  outside  computa¬ 
tional  boundaries--two  cards  must  be 
included 

ZB  Heights  (ft)  for  barrier  cells  (see  card 

group  17)--start  on  a  new  card--two  cards 
must  be  included 

CB  Chezy  coefficients  for  barrier  cells  (see 

card  group  17)--start  on  a  new  card--two 
cards  must  be  included 


12  XMAN. 

(10F8.3)  1 


38 


,m*  +  ***  i 


^  A  ^.mi  .4,4  mA  ■ 


Necessity 

Code 


Card  Group 
(Format) 


Variable 


Description 


13  N  Corrections  to  bed  elevation  matrix — grid 

(215, F8. 3)  M  indices  N,M  and  corrected  bed  elevation 

DNM  DNM  (ft)--NZP  cards  must  be  provided 

Note:  Group  13  is  omitted  if  NZP=0 

14  ISH0RE2^  Matrix  indicating  shoreline  position--for 

(2014)  each  N,  read  an  M  value  corresponding  to 

the  last  (seaward)  land  cell  between 
beach  and  ocean.  There  are  NMAX  values 

IBRKl^j  Matrix  indicating  first  breaker  line 

position--for  each  N,  read  an  M  value  cor¬ 
responding  to  the  breaker  line 

ISH0RE3jj  Matrix  indicating  a  fictitious  shoreline 
position  for  a  second  surf  zone 

111111(2^  Matrix  indicating  second  breaker  line 
position 

Note:  If  there  is  no  second  surf  zone,  set  all  the 

values  of  1SH0RE3  and  IBRK2  to  zero 


15 

(3512) 


17 

(312,414) 


Friction  codes  (1  to  20) — for  each  N, 
read  MANN  M,  M  =  1,  MMAX  (start  a  new 
card  for  £ach  N) — use  special  code  of  99 
for  all  water  cells  outside  of  computa¬ 
tional  boundaries — code  99  corresponds  to 
XMAN(l) 

Used  for  boundary  data 
1  —  impermeable  nonovertopping  barrier 

8 —  elevation  (tidal)  input  boundary 
(used  for  radiation  boundary) 

9- -flow  input  boundary  (used  for  uniform 

flux  boundary) 

99 — end  of  group  17  data 

0--for  flux  boundary 

l,2,3--used  for  indicating  different  ele¬ 
vation  (tidal)  boundaries 

1 —  if  flow  or  elevation  (tide)  is  in  the 
x-direction 

2- -if  flow  or  elevation  (tide)  is  in  the 

y-direction 

Grid  index  (M  or  N  value)  of  boundary 
line 

Boundary  extends  from  N  or  M  =  12  to  13 


•  '  •  '  ■  »  »  »  •  *  *  •  m  •  •  •  •  ^  «  «  »  •  •  »  •  “  »  •  •  *  '  •  • 

•  „  •  .  <  .  *  •  *  *  *  «,*  »  •  “  V  ■*  %  *.*  «  •»*  «.*  S*  ■S.'*  •  '  •  ‘  *,  %  •  *  *.  %  \  '  »,  ‘  % 
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Necessity 

Code 


R 

(Continued) 


0 


R 

R 


Card  Group 

(Format)  Variable  _ Description _ 

17  14  0 — for  top  or  left  boundary  (low  values 

(312,414)  of  x  or  y) 

1 — for  bottom  or  right  boundary  (high 
values  of  x  or  y) 

Note:  The  termination  of  boundary  data  (group  17) 
is  indicated  by  including  at  the  end  of  the 
group  a  card  containing  99  in  columns  1  to  2. 
All  other  columns  should  be  blank 


20  SURIN.  Entries  in  input  elevation  (tidal) 

(10F8.4)  1  tables--there  are  NTID  tables  and  for 

each  table,  there  are  ITID  entries--the 
order  in  which  the  tables  are  read  in  de¬ 
pends  on  the  sequence  of  the  elevation 
(tidal)  boundaries  in  group  17--start 
each  table  on  a  new  card 


Note:  For  radiation  boundary,  set  ITID  =  2,  JTID  >  MAXTIM 
and  SURIN j  =  SURIN2  =0. 

Note:  Group  20  is  omitted  if  NTID  =  0 


23 

JDELAY 

Delay  saving  gage  information  until 

(1615) 

ITIME  =  JDELAY 

32 

RHO 

3 

Density  of  seawater  (slugs/ft  ) — usually 

(NAMELIST 

set  to  1.99 

$PAR) 

EY 

2 

Eddy  viscosity  ey  (ft  /sec) 

RAD 

A  weighting  factor  (between  zero  and  one) 

NTIMEB 

for  radiation  stress  gradient  terms 

Number  of  time-steps  over  which  radiation 

ADV1 

terms  are  to  be  built  up 

A  weighting  factor  (between  zero  and  one) 

FRC1 

for  advection  terms 

A  weighting  factor  (between  zero  and  one) 

ANLH 

for  bottom  friction  terms 

The  parameter  NT„  in  Longuet-Higgins  ex¬ 
pression  for  eday  viscosity  (Equa¬ 
tion  13)--varies  between  0  and  0.016 

Drag  coefficient  c  in  the  Longuet- 

CF 

ITAVG 

Higgins  linear  relation  for  bed  friction 

Number  of  time-steps  over  which  the  solu¬ 

tion  (values  of  n  ,  U  ,  V  ,  and  dis¬ 
charges)  is  to  be  averaged.  The  program 
averages  the  solution  starting  with 
values  from  step  MAXTIM-ITAVG  +  1.  These 
values  are  printed  as  well  as  saved  on 
appropriate  files 

S*. 


It 


40 


Necessity  Card  Group 


_ Description _ 

Index  which  indicates  method  of  computing 
eddy  viscosity  e 
1  =  Longuet-Higgins  method 
1  =  Jonsson  et  al.  method 
3  =  eddy  viscosity  £  is  read  in  from 
file  40  X 

The  matrix  containing  the  distance  from 
the  shoreline  to  each  cell  for  use  in 
Longuet-Higgins  formula  for  eddy  viscosity 
e  .  For  each  value  of  N,  read  MMAX  val- 
ues.  Start  each  new  N  on  a  separate  card. 
Note:  For  water  cells  offshore  of  breaker 
line,  keep  XDIST  same  as  the  value 
at  the  breaker  line 

Input  Files 

42.  In  Boeing  Computer  Services  (BCS),  the  source  deck  for  the  program 
is  on  file  WIFMSRC.  It  has  to  be  updated  first  and  a  program  library  created. 
Then  the  program  library  should  be  updated  again  with  the  update  file  UPDL2RV. 

43.  The  following  input  data  files  are  often  used.  Some  are  optional 
depending  on  the  run. 

File  FT  35 

44.  This  file  contains  the  wave  data,  usually  the  output  from  a  numeri¬ 
cal  wave  climate  program.  For  each  grid  cell,  the  wave  number  AK,  wave  angle 
TH,  and  the  wave  height  HT  are  read  in  sequence.  The  following  statements  are 
used  to  read  the  file  in  subroutine  DEPTH: 

DO  10060  M=1 ,  MMAX 

DO  10060  N=l,  NMAX 

I  =  NMAX*(M- 1 )  +  N 

READ  (35,10061)  AK(I),  TH(I),  HT(I) 

10060  CONTINUE 

10061  FORMAT  (F10.6,  2F10.3) 

Note  that  in  the  program  the  same  matrix  may  be  used  either  as  a  double  index, 
e.g.,  AK(N,M)  or  a  single  index,  e.g.,  AK(I)  array.  The  conversion  from  the 
two  index  to  one  index  system  is  accomplished  by  the  relation 
I  =  NMAX*(M-1)  ♦  N 


Code 

(Format) 

Variable 

R 

33 

I  EDDY 

(12) 

0 

34 

(16F5.0) 

®ISTN.M 

45.  This  file  contains  the  bed  elevations  D(N,M)  for  grid  cells.  The 
following  statements  are  used  to  read  the  file  in  subroutine  DEPTH: 

DO  99005  M=1 ,MMAX 

READ  (36,99001)  (D(N,M) ,  N  =  l.NMAX) 

99001  FORMAT  (10F8.3) 

99005  CONTINUE 

Note:  Usually,  if  SWL  is  taken  as  the  datum, 

D(N,M)  <  0  for  water  cells 
>  0  for  land  cells 
This  is  a  required  file. 

File  FT  39 

46.  This  file  contains  the  grid  expansion  coefficients  YNU  and  XMU  for 
a  variable  grid.  They  correspond  to  (J^  an<*  Mj  »  respectively.  There  will 
be  NYY  =  2*  NMAX  values  of  YNU  and  NXX  =  2*MMAX  values  of  XMU.  These  coef¬ 
ficients  are  provided  by  a  special  program  called  GRID.  The  following  state¬ 
ments  are  used  to  read  the  file  in  subroutine  INTAKE: 

READ  (39,99002)  (YNU(I),  I  =  1 ,NYY) 

READ  (39,99002)  (XMU(I),  I  =  1 ,NXX) 

99002  FORMAT  (4G20.11) 

This  file  is  optional  and  required  only  for  a  variable  grid  (IXPAN  f  0). 

File  FT  40 

47.  This  file  contains  the  eddy  viscosity  coefficients  for  the 
x-direction,  EX.  The  following  statements  are  used  to  read  the  coefficients 
in  a  string  in  subroutine  EDDYVIS: 

READ  (40,115)  (EX(I),I  =  1,NMX) 

115  FORMAT  (16F5.1) 

This  file  is  optional  and  required  only  if  IEDDY  =  3. 

File  FT  95 

48.  This  is  a  "card  image"  file.  It  contains  the  rest  of  the  input 
data  needed  for  running  the  program.  The  number  of  the  file  (NDTAP)  may  be 
other  than  95.  The  number  should  be  indicated  always  on  the  first  input  data 
card  (card  group  1A) .  This  file  is  required . 


PART  VI:  MODEL  OUTPUT 


General  Description 

49.  This  section  describes  the  general  output  in  the  form  of  printed 
output  and  output  files  that  may  be  obtained  from  the  model  and  the  controls 
one  may  exercise  on  the  same.  Most  of  the  input  data  are  printed  as  soon  as 
they  are  read  in  so  that  errors  in  the  input  data  may  be  easily  detected  and 
corrected  by  the  user.  By  setting  the  variable  MPR  (card  group  3)  appropri¬ 
ately,  the  user  may  obtain  as  much  information  as  he  or  she  desires  on  the  flag 
arrays,  depths,  flood,  barrier,  and  tidal  data.  Similarly,  by  using  the  NPRINT 
(card  group  7)  option,  the  user  may  obtain  detailed  printout  of  q  ,  U  ,  and 

V  over  the  whole  grid  at  intermediate  times  during  the  computation.  By  using 
the  variables  NGAGE,  NFREQ,  NPOT,  and  MPOT  (card  groups  5  and  8),  the  user  may 
provide  for  special  gages  at  several  locations  within  the  grid  and  obtain  a 
time-history  of  water  levels  and  velocities  there  either  to  check  the  results 
of  the  model  or  to  compare  model  results  with  actual  field  gage  data.  The  user 
must  provide  the  necessary  Job  Control  Language  (JCL)  to  print  the  files  FT  07 
and  FT  08.  By  setting  the  variable  ITAVG  (card  group  32),  the  user  may  obtain 
a  solution  that  is  averaged  over  the  last  ITAVG  time-steps.  The  program  prints 
the  averaged  water  levels  and  velocity  components  at  the  cell  centers  as  well 
as  the  discharges  across  cell  faces  and  stores  them  on  file  FT  20  for  later 
use  in  a  sediment  transport  model. 

50.  Program  CURRENT  constructs  various  arrays  and  tables  packed  with 

data  to  control  double-sweep  computation,  forcing  boundaries,  barrier  treat¬ 
ment,  etc.  There  are  two  flag  arrays  ICU(N,M)  and  ICV(N,M)  to  control  computa¬ 
tion  in  x-  and  y-directions ,  respectively.  Each  element  of  ICU  (ICV)  consists 
of  two  digits:  njn2  •  When  the  flags  are  printed,  for  any  given  cell  (N,M) 
the  values  of  ICU  and  ICV  are  printed  together  as  a  four-digit  number  ICF.  The 
first  digit  in  ICU  (ICV)  defines  the  character  of  the  cell  on  the  bottom 

(right)  face  of  the  cell.  Digits  n^  and  n^  have  the  following  significance; 

nl  _ Definition _ 

1  Bottom  (right)  face  of  cell  is  an  exposed  barrier  that  cannot  be  over¬ 
topped.  n^  has  no  meaning.  Note:  The  code  ICU=10,  ICV=10  is  set  by 
the  program  to  indicate  water  cells  outside  computational  boundaries. 

(Continued) 


Definition 


5  No  flow  through  bottom  (right)  face  of  cell,  n^  is  set  to  zero. 

6  Flow  through  bottom  (right)  face  of  cell,  n ^  indicates  the  type  of 
computation  to  use  for  advection.  =  0  means  no  advection,  n^  =  1 
advection  in  x-direction  only,  n^  =  2  advection  in  y-direction  only,  and 
n2  =  3  advection  in  both  x-  and  y-directions . 

7  Constructed  by  the  program  to  indicate  no  computation  for  x-(y-)  sweep. 
n2  is  set  to  zero  to  indicate  a  lower  boundary,  and  1  for  an  upper 
boundary,  for  the  other  sweep  direction. 

8  Elevation  (tidal)  or  radiation  boundary  cell.  n2  indicates  the  number 
of  the  forcing  elevation  (tide)  used. 

9  Bottom  (right)  face  is  a  discharge  boundary,  =  0  indicates  a  uniform 

flux  boundary  condition  applied  at  cell  center;  otherwise  indicates 

the  number  of  the  forcing  discharge  used. 


Forcing  elevation  (tidal)  boundary  control  vectors  take  the  form: 

Vector  i  :  1000000*INDX+1000*N+M 

and  discharge  boundary  vectors  the  form: 

Vector  i  :  1000000*(INDX+10*IDIR)+1000*N+M 

where  i  is  the  ith  element  of  the  vector,  N  and  M  are  indices  of  the  grid 
cell  and  INDX  and  IDIR  are  defined  in  card  group  17.  The  user  is  urged  to  set 
MAXTIM=0  (card  group  5)  and  MPR=-1  (card  group  3)  during  the  first  run  to 
check  the  accuracy  of  the  input  data,  the  flags,  boundary  vectors,  etc.  Once 
these  are  found  to  be  correct,  then  MAXTIM  and  MPR  may  be  set  to  desired 
values  and  the  actual  computations  may  be  performed. 


Output  Files 

51.  Output  may  be  stored  by  the  model  in  the  following  files.  Some  are 
optional.  The  user  must  dispose  the  files  with  proper  JCL  if  he  or  she  wants 
to  save  the  files  for  later  use. 

File  FT  01 

52.  This  file  contains  the  bed  elevations  h  ,  and  rj  ,  U  ,  and  V 
saved  every  INTAP  time-steps  in  subroutine  POUT.  The  reader  is  referred  to 
that  subroutine  for  details. 

Files  FT  07  and  08 

53.  These  contain  the  results  n  ,  U  ,  and  V  for  selected  gages  (see 


card  groups  5  and  8)  stored  at  intervals  of  NFREQ ,  starting  from  time-step 
JDELAY.  Each  file  contains  information  for  up  to  10  gages.  It  is  convenient 
to  copy  these  files  to  system  OUTPUT  file  at  the  end  of  the  run.  The  results 
will  then  be  printed  in  a  tabular  form.  Files  FT  07  and  FT  08  are  written  on 
in  subroutine  DRIVE. 

File  FT  11 

54.  This  file  contains  the  surface  elevation  rj  and  the  velocities 
U  ,  V  (averaged,  respectively,  at  the  center  of  the  cell)  for  the  special 
gage  locations  in  a  form  that  is  convenient  for  plotting  time  series.  The 
results  are  stored  at  intervals  of  NFREQ,  starting  from  time-step  JDELAY, 
using  the  following  statements: 

DO  10020  J=l,  NGAGE 
N  =  NPOT(J) 

M  =  MPOT(J) 

II  =  NMAX*(M-1)+N 

WRITE(11)ITIME,J,N,M,SEP(II),UU,W 
10020  CONTINUE 

where  ITIME  is  the  time-step  number.  This  file  is  written  on  in  subroutine 
DRIVE. 

File  FT  13 

55.  This  file  contains  the  surface  elevations  all  over  the  grid  at 
selected  instants  of  time  NPRINT  during  the  course  of  the  computation.  The 
file  is  written  on  in  subroutine  POUT  using  the  following  statements: 

WRITE(13)  ITIME 
DO  4800  M=1  ,MMAX 

DO  4800  N=1 ,NMAX 

II  =  (M-1)*NMAX+N 
WRITE (13)SEP(I1) 

4800  CONTINUE 

This  file  may  be  used  for  3-D  plotting. 

File  FT  14 

56.  This  file  contains  the  velocities  U  ,  V  at  the  center  of  each 
grid  cell  for  all  the  grid  cells  at  selected  instants  of  time  NPRINT  during 
the  course  of  the  computation.  The  file  is  written  on  in  subroutine  POUT 
using  the  following  statements: 


WRITE ( 14) ITIME 
DO  4900  M=1  ,MMAX 

DO  4900  N=1,NMAX 
II=NMAX*(M-1)+N 
WRITE(14)  UU,W 
4900  CONTINUE 

This  file  may  be  used  for  vector  plotting. 

File  FT  15 

57.  This  file  contains  the  vertically  integrated  discharges  DISCHX  and 
DISCHY  across  the  lower  and  right  faces  of  each  grid  call  for  all  the  grid 
cells  at  the  end  of  the  run  (ITIME=MAXTIM) .  This  file  is  written  on  in  sub¬ 
routine  POUT  using  the  following  statements: 

DO  31000  M=l,  MMAX 

DO  31000  N=1 ,  NMAX 

II  =  NMAX*(M-1)+N 

WRITE ( 15 )DISCHX ( I I ) , DISCHY ( I I ) 

31000  CONTINUE 

File  FT  20 

58.  This  file  contains  some  of  the  information  needed  for  running  a 
numerical  sediment  transport  model  subsequently,  such  as  shoreline  and  breaker 
line  positions,  updated  wave  information,  average  values  (averaged  over  ITAVG 
time-steps)  of  r)  >  U  ,  V  ,  and  discharges  across  cell  faces.  It  is  written 
on  in  subroutines  DEPTH,  M0TN4,  and  POUT,  respectively.  The  reader  is  referred 
to  these  subroutines  for  details. 


PART  VII:  MODEL  APPLICATION 


Plane  Beach:  Longshore  Currents 

59.  In  order  to  demonstrate  the  use  of  the  model  for  a  simple  situa¬ 
tion,  the  case  of  monochromatic  waves  obliquely  incident  on  a  plane  beach  will 
be  considered.  This  is  the  same  case  that  was  discussed  in  paragraphs  27-30 
of  PART  IV  except  that  the  effect  of  setup,  mixing,  and  advection  are  taken 
into  account  now  during  the  course  of  the  computation.  The  beach  has  a  con¬ 
stant  slope  of  1:30.  The  wave  characteristics  in  deep  water  are  given  by 

T  =  12  sec,  =  10  ft,  and  0^  =  20  deg.  For  purposes  of  the  present  Simula 
tion,  a  uniform  grid  with  Ax  =  Ay  =  60  ft,  NMAX  =  6  ,  and  MMAX  =  50  is 
chosen.  Also,  let  c  =  0.01  ,  y  =  0.82  ,  N^H  =  0.00783  ,  At  =  5.0  sec,  and 
NTIMEB  =  15.  Uniform  flux  boundary  conditions  are  employed  at  the  lateral 
boundaries  and  a  radiation  boundary  condition  at  the  offshore  boundary.  The 
simulation  is  run  for  105  time-steps. 

60.  The  CRAY  J CL  for  running  the  program  at  BCS  is  shown  in  Figure  14. 
The  user  must  supply  his  or  her  own  user  number  (UN)  and  password  (PW)  and 
filename  XXXLOG  for  dayfile.  File  WAVESRV  has  the  results  of  a  wave  climate 
program  for  the  case  under  consideration.  File  DEPBCH  has  the  bathymetric  in¬ 
formation  and  file  DATABCH  most  of  the  input  data.  No  output  files  are  saved. 

61.  The  input  data  for  the  job  from  file  DATABCH  are  shown  in  Figure  15 
Note  that  the  full  print  option  MPR=-1  has  been  selected  and  NGAGE  has  been 
set  to  20  with  NFREQ=5  and  JDELAY=0.  As  for  eddy  viscosity,  IEDDY  has  been 
set  to  1  (Longuet-Higgins  formulation)  so  that  both  ANLH  and  XDIST  have  to  be 
supplied.  Note  also  that  blank  lines  are  used  for  input  statements  containing 
all  zeros. 

62.  Figure  16  shows  the  "echo-print"  of  the  input  data  by  the  program, 
Figure  17  the  results  of  computation,  and  Figure  18  a  sample  printing  of  the 
results  for  selected  gages.  Note  that  the  program  prints  the  updated  wave 
heights  at  the  end  of  the  run  and  the  CPU  time  spent  in  the  double-sweep 
computation. 

Oregon  Inlet 

63.  The  next  application  to  be  considered  is  that  of  Oregon  Inlet. 

This  has  been  discussed  at  length  in  paragraphs  32-37  of  PART  IV.  This  is  a 


47 


CDEXX *  P  OS * T4 0 *  STCft 1 . 

USER* UN- PM. 
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FETCH*  DN=FT95* GDN=DPTBBCH* UN=CEROD2*  DT=C» DS=FF. 
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♦MEDR 
♦REhB  8 
♦MEDR 
95 
♦WEOF 

Figure  14.  Job  Control  Language  (JCL)  for  plane  beach  application 
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Figure  15.  Input  data  for  plane  beach  application 
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Figure  16.  Echo  print  of  input  data  for  plane  beach  application 
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Figure  18.  Sample  results  for  selected  gages  for  plane  beach  computation 


more  complicated  situation.  A  variable  grid  is  employed,  with  NMAX=77  and 
MMAX=54,  so  as  to  have  a  greater  resolution  near  the  inlet  and  the  surf  zone. 
The  JCL  for  this  case  is  shown  in  Figure  19.  In  addition  to  files  FT  35,  36, 
and  95,  file  FT  39  containing  the  grid  expansion  coefficients  is  needed.  In 
this  case,  output  files  FT  11,  13,  14,  15,  and  20  are  being  saved  at  the  end 
of  the  run  for  later  use. 

64.  Part  of  the  input  data  from  file  DATAW3L  is  shown  in  Figure  20. 

Note  that  MPR  is  set  to  -1  for  full  printout,  Aor^  =  =  100  ft  and  At 

=  18  sec.  The  simulation  is  run  for  67  time-steps.  Twenty  special  gages  are 
used  with  NFREQ=2.  A  printout  of  results  is  desired  at  ITIME=67.  There  are 
two  surf  zones  over  a  part  of  the  grid.  JDELAY  is  set  to  15  to  delay  printing 
of  gage  information.  It  is  desired  to  build  up  the  radiation  stress  terms 
over  15  time-steps  (NTIMEB=15).  y  is  set  t0  0.754.  Averaging  of  results 
over  the  last  10  time-steps  is  required  (ITAVG=10).  Jonsson-type  eddy  vis¬ 
cosity  formulation  is  desired  so  IEDDY  is  set  to  2. 

65.  Figure  21  shows  the  coded  friction  array  and  the  flags  for  this  ap¬ 
plication.  A  sampling  of  the  averaged  results  for  r|  ,  U  ,  and  V  at  the 
end  is  also  shown. 


Cost  of  Computation 

66.  In  general,  the  cost  of  running  the  model  depends  on  the  number  of 
grid  points  used,  the  number  of  time-steps  of  simulation  run,  amount  of  input 
and  output,  the  job  priority,  the  computer,  and  the  computer  vendor  used.  The 
following  total  cost  information  was  for  running  the  model  at  standard  priority 
on  the  CRAY  computer  of  BCS.  For  the  plane  beach  application,  described  in 
paragraphs  59-62,  involving  300  grid  points  and  105  time-steps,  the  total  cost 
was  approximately  $7.  For  the  Oregon  Inlet  application,  described  in  para¬ 
graphs  63-65,  involving  4,158  grid  points  and  67  time-steps,  the  total  cost 
was  approximately  $23.  Therefore  the  computation  costs  for  the  model  may  be 
considered  modest. 
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Figure  19.  Job  Control  Language  (JCL)  for  Oregon  Inlet, 
North  Carolina,  application 
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Figure  21.  Coded  friction  arr^y,  flags,  and  a  sample  of  averaged  results  for 
Oregon  Inlet,  North  Carolina,  application  (Sheet  1  of  26) 
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PART  VIII:  SUMMARY  AND  CONCLUSIONS 


67.  This  report  describes  the  development  of  a  generalized  numerical 
model  for  coastal  currents  induced  by  breaking  waves.  The  model  employs  the 
radiation  stress  approach  of  Longuet-Higgins  and  solves  the  vertically  inte¬ 
grated  equations  of  momentum  and  continuity  using  an  alternating  direction  im¬ 
plicit  scheme.  It  includes  mixing  and  advection  terms. 

68.  The  model  has  the  following  desired  features  in  terms  of  applica¬ 
tion  to  engineering  projects:  the  ability  to  handle  real-life  bathymetries, 
flexibility  in  terms  of  formulation  of  mixing  and  friction  terms,  computa¬ 
tional  efficiency,  and  economy. 

69.  The  model  was  applied  to  a  case  of  normally  incident  waves  on  a 
plane  beach.  The  results  for  setup  matched  the  experimental  data  of  Bowen, 
Inman,  and  Simmons  (1968). 

70.  For  oblique  incidence  of  waves  on  a  plane  beach,  model  results  for 
longshore  currents  were  compared  with  the  analytical  solution  of  Longuet- 
Higgins,  first  neglecting  the  effect  of  setup  and  later  including  the  effect 
of  setup.  Agreement  was  excellent.  As  the  numerical  grid  was  made  finer,  the 
numerical  results  tended  to  converge  toward  the  analytical  solution.  Next, 
the  effect  of  lateral  mixing  on  the  velocity  distribution  was  studied.  As  the 
mixing  parameter  P  increased,  the  numerical  solution  exhibited  the  proper 
behavior,  since  the  magnitude  of  the  peak  decreased,  the  peak  moved  closer  to 
the  shoreline,  and  the  velocities  offshore  of  the  breaker  line  increased. 

71.  The  model  was  finally  applied  to  a  field  situation  corresponding  to 
Oregon  Inlet,  North  Carolina.  The  bathymetry  was  very  irregular  and  complex 
owing  to  the  presence  of  channels  and  shoals.  A  variable  grid  was  used.  The 
significant  wave  condition  during  a  part  of  the  Ash  Wednesday  storm  of  March 
1962  was  selected  for  simulation.  The  numerical  results  obtained  for  this 
case  appeared  to  be  reasonable  and  the  computer  costs  were  modest. 

72.  For  the  convenience  of  the  potential  user,  model  input,  output,  and 
files  are  described  and  two  sample  applications  are  presented. 
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APPENDIX  A:  NOTATION 
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Mapping  function,  dimensionless 
Coefficients,  1/sec 
Coefficients,  sec 
Coefficient,  dimensionless 
Computational  parameter,  ft 
Mapping  function,  dimensionless 
Computational  parameter,  ft/sec 

Drag  coefficient  (of  the  order  of  0.01),  dimensionless 
Mapping  function,  dimensionless 
Phase  speed,  ft/sec 

Total  water  depth  including  setup,  ft 
Water  depth  at  wave  breaking,  ft 
Wave  energy  density,  lb/ft 

2 

Acceleration  due  to  gravity,  32.174  ft/sec 

Bed  elevation  with  still-water  level  taken  as  zero,  ft 

Local  still-water  depth,  ft 

Local  wave  height,  ft 

Breaking  wave  height,  ft 

Deepwater  wave  height,  ft 

Local  wave  number,  2n/L  ,  1/ft 

Local  wavelength,  ft 

Breaking  wavelength,  ft 

x-index  of  arbitrary  cell  center 

y-index  of  arbitrary  cell  center 

Ratio  of  wave  group  celerity  to  wave  phase  celerity, 
dimensionless 

Empirical  coefficient  that  varies  from  0.000  to  0.016, 
dimensionless 

Mixing  parameter  which  varies  between  0  and  0.40, 
dimensionless 

Recursion  coefficient,  sec 
Recursion  coefficient,  sec 
Recursion  coefficient,  ft 
Recursion  coefficient,  ft 


Recursion  coefficient,  1/sec 
Recursion  coefficient,  1/sec 
Beach  slope,  dimensionless 
Recursion  coefficient,  ft/sec 
Recursion  coefficient,  ft/sec 

Radiation  stress  in  the  x-direction  (normal  to  the 
y-z  plane),  lb/ft 

Radiation  stress  in  the  y-direction  in  the  x-z  plane,  lb/ft 

Radiation  stress  in  the  y-direction  (normal  to  the 
x-z  plane),  lb/ft 

Time,  sec 

Wave  period,  sec 

Computational  parameter,  dimensionless 
Computational  parameter,  dimensionless 

Depth-averaged  horizontal  velocity  in  x-direction,  ft/sec 

Time-averaged  wave  orbital  velocity  at  the  bottom,  ft/sec 

Depth-averaged  horizontal  velocity  in  y-direction,  ft/sec 

Horizontal  Cartesian  coordinate,  ft 

Horizontal  Cartesian  coordinate,  ft 

Vertical  Cartesian  coordinate,  ft 

Arbitrary  variable 

Computational  space  coordinate,  ft 

Computational  space  coordinate,  ft 

Wave  breaking  index  on  the  order  of  1.0,  dimensionesss 
Difference  operator 

Time  increment,  sec 

Cell  dimension  in  x-direction  in  real  space,  ft 


Cell  dimension  in  y-direction  in  real  space,  ft 

Cell  dimension  in  o^-direction  in  computational  space,  ft 

Cell  dimension  in  a. -direction  in  computational  space,  ft 

2  2 

Eddy  viscosity  in  the  x-direction,  ft  /sec 

2 

Eddy  viscosity  in  the  x-direction,  ft  /sec 

Water-surface  elevation  at  intermediate  time  level,  ft 

Displacement  of  the  mean  free  surface  with  respect  to  still- 
water  level,  ft 

Local  wave  direction,  deg 
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Symbol 


**1 
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T 

xy 


Expansion  coefficient,  dimensionless 
Expansion  coefficient,  dimensionless 
3.14159,  dimensionless 

2  4 

Mass  density  of  seawater,  1.99  lb-sec  /ft 

2 

Bottom  friction  stress  in  the  x-direction,  lb/ft 

2 

Bottom  friction  stress  in  the  y-direction,  lb/ft 

2 

Lateral  shear  stress  due  to  turbulent  mixing,  lb/ft 


3  Partial  derivative  symbol,  dimensionless 
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